// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_SUPERLUSUPPORT_H
#define EIGEN_SUPERLUSUPPORT_H

namespace Eigen {

#if defined(SUPERLU_MAJOR_VERSION) && (SUPERLU_MAJOR_VERSION >= 5)
#define DECL_GSSVX(PREFIX, FLOATTYPE, KEYTYPE)                                                                                                              \
    extern "C" {                                                                                                                                            \
    extern void PREFIX##gssvx(superlu_options_t*,                                                                                                           \
                              SuperMatrix*,                                                                                                                 \
                              int*,                                                                                                                         \
                              int*,                                                                                                                         \
                              int*,                                                                                                                         \
                              char*,                                                                                                                        \
                              FLOATTYPE*,                                                                                                                   \
                              FLOATTYPE*,                                                                                                                   \
                              SuperMatrix*,                                                                                                                 \
                              SuperMatrix*,                                                                                                                 \
                              void*,                                                                                                                        \
                              int,                                                                                                                          \
                              SuperMatrix*,                                                                                                                 \
                              SuperMatrix*,                                                                                                                 \
                              FLOATTYPE*,                                                                                                                   \
                              FLOATTYPE*,                                                                                                                   \
                              FLOATTYPE*,                                                                                                                   \
                              FLOATTYPE*,                                                                                                                   \
                              GlobalLU_t*,                                                                                                                  \
                              mem_usage_t*,                                                                                                                 \
                              SuperLUStat_t*,                                                                                                               \
                              int*);                                                                                                                        \
    }                                                                                                                                                       \
    inline float SuperLU_gssvx(superlu_options_t* options,                                                                                                  \
                               SuperMatrix* A,                                                                                                              \
                               int* perm_c,                                                                                                                 \
                               int* perm_r,                                                                                                                 \
                               int* etree,                                                                                                                  \
                               char* equed,                                                                                                                 \
                               FLOATTYPE* R,                                                                                                                \
                               FLOATTYPE* C,                                                                                                                \
                               SuperMatrix* L,                                                                                                              \
                               SuperMatrix* U,                                                                                                              \
                               void* work,                                                                                                                  \
                               int lwork,                                                                                                                   \
                               SuperMatrix* B,                                                                                                              \
                               SuperMatrix* X,                                                                                                              \
                               FLOATTYPE* recip_pivot_growth,                                                                                               \
                               FLOATTYPE* rcond,                                                                                                            \
                               FLOATTYPE* ferr,                                                                                                             \
                               FLOATTYPE* berr,                                                                                                             \
                               SuperLUStat_t* stats,                                                                                                        \
                               int* info,                                                                                                                   \
                               KEYTYPE)                                                                                                                     \
    {                                                                                                                                                       \
        mem_usage_t mem_usage;                                                                                                                              \
        GlobalLU_t gLU;                                                                                                                                     \
        PREFIX##gssvx(                                                                                                                                      \
            options, A, perm_c, perm_r, etree, equed, R, C, L, U, work, lwork, B, X, recip_pivot_growth, rcond, ferr, berr, &gLU, &mem_usage, stats, info); \
        return mem_usage.for_lu; /* bytes used by the factor storage */                                                                                     \
    }
#else  // version < 5.0
#define DECL_GSSVX(PREFIX, FLOATTYPE, KEYTYPE)                                                                                                        \
    extern "C" {                                                                                                                                      \
    extern void PREFIX##gssvx(superlu_options_t*,                                                                                                     \
                              SuperMatrix*,                                                                                                           \
                              int*,                                                                                                                   \
                              int*,                                                                                                                   \
                              int*,                                                                                                                   \
                              char*,                                                                                                                  \
                              FLOATTYPE*,                                                                                                             \
                              FLOATTYPE*,                                                                                                             \
                              SuperMatrix*,                                                                                                           \
                              SuperMatrix*,                                                                                                           \
                              void*,                                                                                                                  \
                              int,                                                                                                                    \
                              SuperMatrix*,                                                                                                           \
                              SuperMatrix*,                                                                                                           \
                              FLOATTYPE*,                                                                                                             \
                              FLOATTYPE*,                                                                                                             \
                              FLOATTYPE*,                                                                                                             \
                              FLOATTYPE*,                                                                                                             \
                              mem_usage_t*,                                                                                                           \
                              SuperLUStat_t*,                                                                                                         \
                              int*);                                                                                                                  \
    }                                                                                                                                                 \
    inline float SuperLU_gssvx(superlu_options_t* options,                                                                                            \
                               SuperMatrix* A,                                                                                                        \
                               int* perm_c,                                                                                                           \
                               int* perm_r,                                                                                                           \
                               int* etree,                                                                                                            \
                               char* equed,                                                                                                           \
                               FLOATTYPE* R,                                                                                                          \
                               FLOATTYPE* C,                                                                                                          \
                               SuperMatrix* L,                                                                                                        \
                               SuperMatrix* U,                                                                                                        \
                               void* work,                                                                                                            \
                               int lwork,                                                                                                             \
                               SuperMatrix* B,                                                                                                        \
                               SuperMatrix* X,                                                                                                        \
                               FLOATTYPE* recip_pivot_growth,                                                                                         \
                               FLOATTYPE* rcond,                                                                                                      \
                               FLOATTYPE* ferr,                                                                                                       \
                               FLOATTYPE* berr,                                                                                                       \
                               SuperLUStat_t* stats,                                                                                                  \
                               int* info,                                                                                                             \
                               KEYTYPE)                                                                                                               \
    {                                                                                                                                                 \
        mem_usage_t mem_usage;                                                                                                                        \
        PREFIX##gssvx(                                                                                                                                \
            options, A, perm_c, perm_r, etree, equed, R, C, L, U, work, lwork, B, X, recip_pivot_growth, rcond, ferr, berr, &mem_usage, stats, info); \
        return mem_usage.for_lu; /* bytes used by the factor storage */                                                                               \
    }
#endif

DECL_GSSVX(s, float, float)
DECL_GSSVX(c, float, std::complex<float>)
DECL_GSSVX(d, double, double)
DECL_GSSVX(z, double, std::complex<double>)

#ifdef MILU_ALPHA
#define EIGEN_SUPERLU_HAS_ILU
#endif

#ifdef EIGEN_SUPERLU_HAS_ILU

// similarly for the incomplete factorization using gsisx
#define DECL_GSISX(PREFIX, FLOATTYPE, KEYTYPE)                                                                                                      \
    extern "C" {                                                                                                                                    \
    extern void PREFIX##gsisx(superlu_options_t*,                                                                                                   \
                              SuperMatrix*,                                                                                                         \
                              int*,                                                                                                                 \
                              int*,                                                                                                                 \
                              int*,                                                                                                                 \
                              char*,                                                                                                                \
                              FLOATTYPE*,                                                                                                           \
                              FLOATTYPE*,                                                                                                           \
                              SuperMatrix*,                                                                                                         \
                              SuperMatrix*,                                                                                                         \
                              void*,                                                                                                                \
                              int,                                                                                                                  \
                              SuperMatrix*,                                                                                                         \
                              SuperMatrix*,                                                                                                         \
                              FLOATTYPE*,                                                                                                           \
                              FLOATTYPE*,                                                                                                           \
                              mem_usage_t*,                                                                                                         \
                              SuperLUStat_t*,                                                                                                       \
                              int*);                                                                                                                \
    }                                                                                                                                               \
    inline float SuperLU_gsisx(superlu_options_t* options,                                                                                          \
                               SuperMatrix* A,                                                                                                      \
                               int* perm_c,                                                                                                         \
                               int* perm_r,                                                                                                         \
                               int* etree,                                                                                                          \
                               char* equed,                                                                                                         \
                               FLOATTYPE* R,                                                                                                        \
                               FLOATTYPE* C,                                                                                                        \
                               SuperMatrix* L,                                                                                                      \
                               SuperMatrix* U,                                                                                                      \
                               void* work,                                                                                                          \
                               int lwork,                                                                                                           \
                               SuperMatrix* B,                                                                                                      \
                               SuperMatrix* X,                                                                                                      \
                               FLOATTYPE* recip_pivot_growth,                                                                                       \
                               FLOATTYPE* rcond,                                                                                                    \
                               SuperLUStat_t* stats,                                                                                                \
                               int* info,                                                                                                           \
                               KEYTYPE)                                                                                                             \
    {                                                                                                                                               \
        mem_usage_t mem_usage;                                                                                                                      \
        PREFIX##gsisx(options, A, perm_c, perm_r, etree, equed, R, C, L, U, work, lwork, B, X, recip_pivot_growth, rcond, &mem_usage, stats, info); \
        return mem_usage.for_lu; /* bytes used by the factor storage */                                                                             \
    }

DECL_GSISX(s, float, float)
DECL_GSISX(c, float, std::complex<float>)
DECL_GSISX(d, double, double)
DECL_GSISX(z, double, std::complex<double>)

#endif

template <typename MatrixType> struct SluMatrixMapHelper;

/** \internal
  *
  * A wrapper class for SuperLU matrices. It supports only compressed sparse matrices
  * and dense matrices. Supernodal and other fancy format are not supported by this wrapper.
  *
  * This wrapper class mainly aims to avoids the need of dynamic allocation of the storage structure.
  */
struct SluMatrix : SuperMatrix
{
    SluMatrix() { Store = &storage; }

    SluMatrix(const SluMatrix& other) : SuperMatrix(other)
    {
        Store = &storage;
        storage = other.storage;
    }

    SluMatrix& operator=(const SluMatrix& other)
    {
        SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
        Store = &storage;
        storage = other.storage;
        return *this;
    }

    struct
    {
        union
        {
            int nnz;
            int lda;
        };
        void* values;
        int* innerInd;
        int* outerInd;
    } storage;

    void setStorageType(Stype_t t)
    {
        Stype = t;
        if (t == SLU_NC || t == SLU_NR || t == SLU_DN)
            Store = &storage;
        else
        {
            eigen_assert(false && "storage type not supported");
            Store = 0;
        }
    }

    template <typename Scalar> void setScalarType()
    {
        if (internal::is_same<Scalar, float>::value)
            Dtype = SLU_S;
        else if (internal::is_same<Scalar, double>::value)
            Dtype = SLU_D;
        else if (internal::is_same<Scalar, std::complex<float>>::value)
            Dtype = SLU_C;
        else if (internal::is_same<Scalar, std::complex<double>>::value)
            Dtype = SLU_Z;
        else
        {
            eigen_assert(false && "Scalar type not supported by SuperLU");
        }
    }

    template <typename MatrixType> static SluMatrix Map(MatrixBase<MatrixType>& _mat)
    {
        MatrixType& mat(_mat.derived());
        eigen_assert(((MatrixType::Flags & RowMajorBit) != RowMajorBit) && "row-major dense matrices are not supported by SuperLU");
        SluMatrix res;
        res.setStorageType(SLU_DN);
        res.setScalarType<typename MatrixType::Scalar>();
        res.Mtype = SLU_GE;

        res.nrow = internal::convert_index<int>(mat.rows());
        res.ncol = internal::convert_index<int>(mat.cols());

        res.storage.lda = internal::convert_index<int>(MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride());
        res.storage.values = (void*)(mat.data());
        return res;
    }

    template <typename MatrixType> static SluMatrix Map(SparseMatrixBase<MatrixType>& a_mat)
    {
        MatrixType& mat(a_mat.derived());
        SluMatrix res;
        if ((MatrixType::Flags & RowMajorBit) == RowMajorBit)
        {
            res.setStorageType(SLU_NR);
            res.nrow = internal::convert_index<int>(mat.cols());
            res.ncol = internal::convert_index<int>(mat.rows());
        }
        else
        {
            res.setStorageType(SLU_NC);
            res.nrow = internal::convert_index<int>(mat.rows());
            res.ncol = internal::convert_index<int>(mat.cols());
        }

        res.Mtype = SLU_GE;

        res.storage.nnz = internal::convert_index<int>(mat.nonZeros());
        res.storage.values = mat.valuePtr();
        res.storage.innerInd = mat.innerIndexPtr();
        res.storage.outerInd = mat.outerIndexPtr();

        res.setScalarType<typename MatrixType::Scalar>();

        // FIXME the following is not very accurate
        if (int(MatrixType::Flags) & int(Upper))
            res.Mtype = SLU_TRU;
        if (int(MatrixType::Flags) & int(Lower))
            res.Mtype = SLU_TRL;

        eigen_assert(((int(MatrixType::Flags) & int(SelfAdjoint)) == 0) && "SelfAdjoint matrix shape not supported by SuperLU");

        return res;
    }
};

template <typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols> struct SluMatrixMapHelper<Matrix<Scalar, Rows, Cols, Options, MRows, MCols>>
{
    typedef Matrix<Scalar, Rows, Cols, Options, MRows, MCols> MatrixType;
    static void run(MatrixType& mat, SluMatrix& res)
    {
        eigen_assert(((Options & RowMajor) != RowMajor) && "row-major dense matrices is not supported by SuperLU");
        res.setStorageType(SLU_DN);
        res.setScalarType<Scalar>();
        res.Mtype = SLU_GE;

        res.nrow = mat.rows();
        res.ncol = mat.cols();

        res.storage.lda = mat.outerStride();
        res.storage.values = mat.data();
    }
};

template <typename Derived> struct SluMatrixMapHelper<SparseMatrixBase<Derived>>
{
    typedef Derived MatrixType;
    static void run(MatrixType& mat, SluMatrix& res)
    {
        if ((MatrixType::Flags & RowMajorBit) == RowMajorBit)
        {
            res.setStorageType(SLU_NR);
            res.nrow = mat.cols();
            res.ncol = mat.rows();
        }
        else
        {
            res.setStorageType(SLU_NC);
            res.nrow = mat.rows();
            res.ncol = mat.cols();
        }

        res.Mtype = SLU_GE;

        res.storage.nnz = mat.nonZeros();
        res.storage.values = mat.valuePtr();
        res.storage.innerInd = mat.innerIndexPtr();
        res.storage.outerInd = mat.outerIndexPtr();

        res.setScalarType<typename MatrixType::Scalar>();

        // FIXME the following is not very accurate
        if (MatrixType::Flags & Upper)
            res.Mtype = SLU_TRU;
        if (MatrixType::Flags & Lower)
            res.Mtype = SLU_TRL;

        eigen_assert(((MatrixType::Flags & SelfAdjoint) == 0) && "SelfAdjoint matrix shape not supported by SuperLU");
    }
};

namespace internal {

    template <typename MatrixType> SluMatrix asSluMatrix(MatrixType& mat) { return SluMatrix::Map(mat); }

    /** View a Super LU matrix as an Eigen expression */
    template <typename Scalar, int Flags, typename Index> MappedSparseMatrix<Scalar, Flags, Index> map_superlu(SluMatrix& sluMat)
    {
        eigen_assert(((Flags & RowMajor) == RowMajor && sluMat.Stype == SLU_NR) || ((Flags & ColMajor) == ColMajor && sluMat.Stype == SLU_NC));

        Index outerSize = (Flags & RowMajor) == RowMajor ? sluMat.ncol : sluMat.nrow;

        return MappedSparseMatrix<Scalar, Flags, Index>(sluMat.nrow,
                                                        sluMat.ncol,
                                                        sluMat.storage.outerInd[outerSize],
                                                        sluMat.storage.outerInd,
                                                        sluMat.storage.innerInd,
                                                        reinterpret_cast<Scalar*>(sluMat.storage.values));
    }

}  // end namespace internal

/** \ingroup SuperLUSupport_Module
  * \class SuperLUBase
  * \brief The base class for the direct and incomplete LU factorization of SuperLU
  */
template <typename _MatrixType, typename Derived> class SuperLUBase : public SparseSolverBase<Derived>
{
protected:
    typedef SparseSolverBase<Derived> Base;
    using Base::derived;
    using Base::m_isInitialized;

public:
    typedef _MatrixType MatrixType;
    typedef typename MatrixType::Scalar Scalar;
    typedef typename MatrixType::RealScalar RealScalar;
    typedef typename MatrixType::StorageIndex StorageIndex;
    typedef Matrix<Scalar, Dynamic, 1> Vector;
    typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
    typedef Map<PermutationMatrix<Dynamic, Dynamic, int>> PermutationMap;
    typedef SparseMatrix<Scalar> LUMatrixType;
    enum
    {
        ColsAtCompileTime = MatrixType::ColsAtCompileTime,
        MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
    };

public:
    SuperLUBase() {}

    ~SuperLUBase() { clearFactors(); }

    inline Index rows() const { return m_matrix.rows(); }
    inline Index cols() const { return m_matrix.cols(); }

    /** \returns a reference to the Super LU option object to configure the  Super LU algorithms. */
    inline superlu_options_t& options() { return m_sluOptions; }

    /** \brief Reports whether previous computation was successful.
      *
      * \returns \c Success if computation was successful,
      *          \c NumericalIssue if the matrix.appears to be negative.
      */
    ComputationInfo info() const
    {
        eigen_assert(m_isInitialized && "Decomposition is not initialized.");
        return m_info;
    }

    /** Computes the sparse Cholesky decomposition of \a matrix */
    void compute(const MatrixType& matrix)
    {
        derived().analyzePattern(matrix);
        derived().factorize(matrix);
    }

    /** Performs a symbolic decomposition on the sparcity of \a matrix.
      *
      * This function is particularly useful when solving for several problems having the same structure.
      * 
      * \sa factorize()
      */
    void analyzePattern(const MatrixType& /*matrix*/)
    {
        m_isInitialized = true;
        m_info = Success;
        m_analysisIsOk = true;
        m_factorizationIsOk = false;
    }

    template <typename Stream> void dumpMemory(Stream& /*s*/) {}

protected:
    void initFactorization(const MatrixType& a)
    {
        set_default_options(&this->m_sluOptions);

        const Index size = a.rows();
        m_matrix = a;

        m_sluA = internal::asSluMatrix(m_matrix);
        clearFactors();

        m_p.resize(size);
        m_q.resize(size);
        m_sluRscale.resize(size);
        m_sluCscale.resize(size);
        m_sluEtree.resize(size);

        // set empty B and X
        m_sluB.setStorageType(SLU_DN);
        m_sluB.setScalarType<Scalar>();
        m_sluB.Mtype = SLU_GE;
        m_sluB.storage.values = 0;
        m_sluB.nrow = 0;
        m_sluB.ncol = 0;
        m_sluB.storage.lda = internal::convert_index<int>(size);
        m_sluX = m_sluB;

        m_extractedDataAreDirty = true;
    }

    void init()
    {
        m_info = InvalidInput;
        m_isInitialized = false;
        m_sluL.Store = 0;
        m_sluU.Store = 0;
    }

    void extractData() const;

    void clearFactors()
    {
        if (m_sluL.Store)
            Destroy_SuperNode_Matrix(&m_sluL);
        if (m_sluU.Store)
            Destroy_CompCol_Matrix(&m_sluU);

        m_sluL.Store = 0;
        m_sluU.Store = 0;

        memset(&m_sluL, 0, sizeof m_sluL);
        memset(&m_sluU, 0, sizeof m_sluU);
    }

    // cached data to reduce reallocation, etc.
    mutable LUMatrixType m_l;
    mutable LUMatrixType m_u;
    mutable IntColVectorType m_p;
    mutable IntRowVectorType m_q;

    mutable LUMatrixType m_matrix;  // copy of the factorized matrix
    mutable SluMatrix m_sluA;
    mutable SuperMatrix m_sluL, m_sluU;
    mutable SluMatrix m_sluB, m_sluX;
    mutable SuperLUStat_t m_sluStat;
    mutable superlu_options_t m_sluOptions;
    mutable std::vector<int> m_sluEtree;
    mutable Matrix<RealScalar, Dynamic, 1> m_sluRscale, m_sluCscale;
    mutable Matrix<RealScalar, Dynamic, 1> m_sluFerr, m_sluBerr;
    mutable char m_sluEqued;

    mutable ComputationInfo m_info;
    int m_factorizationIsOk;
    int m_analysisIsOk;
    mutable bool m_extractedDataAreDirty;

private:
    SuperLUBase(SuperLUBase&) {}
};

/** \ingroup SuperLUSupport_Module
  * \class SuperLU
  * \brief A sparse direct LU factorization and solver based on the SuperLU library
  *
  * This class allows to solve for A.X = B sparse linear problems via a direct LU factorization
  * using the SuperLU library. The sparse matrix A must be squared and invertible. The vectors or matrices
  * X and B can be either dense or sparse.
  *
  * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
  *
  * \warning This class is only for the 4.x versions of SuperLU. The 3.x and 5.x versions are not supported.
  *
  * \implsparsesolverconcept
  *
  * \sa \ref TutorialSparseSolverConcept, class SparseLU
  */
template <typename _MatrixType> class SuperLU : public SuperLUBase<_MatrixType, SuperLU<_MatrixType>>
{
public:
    typedef SuperLUBase<_MatrixType, SuperLU> Base;
    typedef _MatrixType MatrixType;
    typedef typename Base::Scalar Scalar;
    typedef typename Base::RealScalar RealScalar;
    typedef typename Base::StorageIndex StorageIndex;
    typedef typename Base::IntRowVectorType IntRowVectorType;
    typedef typename Base::IntColVectorType IntColVectorType;
    typedef typename Base::PermutationMap PermutationMap;
    typedef typename Base::LUMatrixType LUMatrixType;
    typedef TriangularView<LUMatrixType, Lower | UnitDiag> LMatrixType;
    typedef TriangularView<LUMatrixType, Upper> UMatrixType;

public:
    using Base::_solve_impl;

    SuperLU() : Base() { init(); }

    explicit SuperLU(const MatrixType& matrix) : Base()
    {
        init();
        Base::compute(matrix);
    }

    ~SuperLU() {}

    /** Performs a symbolic decomposition on the sparcity of \a matrix.
      *
      * This function is particularly useful when solving for several problems having the same structure.
      * 
      * \sa factorize()
      */
    void analyzePattern(const MatrixType& matrix)
    {
        m_info = InvalidInput;
        m_isInitialized = false;
        Base::analyzePattern(matrix);
    }

    /** Performs a numeric decomposition of \a matrix
      *
      * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
      *
      * \sa analyzePattern()
      */
    void factorize(const MatrixType& matrix);

    /** \internal */
    template <typename Rhs, typename Dest> void _solve_impl(const MatrixBase<Rhs>& b, MatrixBase<Dest>& dest) const;

    inline const LMatrixType& matrixL() const
    {
        if (m_extractedDataAreDirty)
            this->extractData();
        return m_l;
    }

    inline const UMatrixType& matrixU() const
    {
        if (m_extractedDataAreDirty)
            this->extractData();
        return m_u;
    }

    inline const IntColVectorType& permutationP() const
    {
        if (m_extractedDataAreDirty)
            this->extractData();
        return m_p;
    }

    inline const IntRowVectorType& permutationQ() const
    {
        if (m_extractedDataAreDirty)
            this->extractData();
        return m_q;
    }

    Scalar determinant() const;

protected:
    using Base::m_l;
    using Base::m_matrix;
    using Base::m_p;
    using Base::m_q;
    using Base::m_sluA;
    using Base::m_sluB;
    using Base::m_sluBerr;
    using Base::m_sluCscale;
    using Base::m_sluEqued;
    using Base::m_sluEtree;
    using Base::m_sluFerr;
    using Base::m_sluL;
    using Base::m_sluOptions;
    using Base::m_sluRscale;
    using Base::m_sluStat;
    using Base::m_sluU;
    using Base::m_sluX;
    using Base::m_u;

    using Base::m_analysisIsOk;
    using Base::m_extractedDataAreDirty;
    using Base::m_factorizationIsOk;
    using Base::m_info;
    using Base::m_isInitialized;

    void init()
    {
        Base::init();

        set_default_options(&this->m_sluOptions);
        m_sluOptions.PrintStat = NO;
        m_sluOptions.ConditionNumber = NO;
        m_sluOptions.Trans = NOTRANS;
        m_sluOptions.ColPerm = COLAMD;
    }

private:
    SuperLU(SuperLU&) {}
};

template <typename MatrixType> void SuperLU<MatrixType>::factorize(const MatrixType& a)
{
    eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
    if (!m_analysisIsOk)
    {
        m_info = InvalidInput;
        return;
    }

    this->initFactorization(a);

    m_sluOptions.ColPerm = COLAMD;
    int info = 0;
    RealScalar recip_pivot_growth, rcond;
    RealScalar ferr, berr;

    StatInit(&m_sluStat);
    SuperLU_gssvx(&m_sluOptions,
                  &m_sluA,
                  m_q.data(),
                  m_p.data(),
                  &m_sluEtree[0],
                  &m_sluEqued,
                  &m_sluRscale[0],
                  &m_sluCscale[0],
                  &m_sluL,
                  &m_sluU,
                  NULL,
                  0,
                  &m_sluB,
                  &m_sluX,
                  &recip_pivot_growth,
                  &rcond,
                  &ferr,
                  &berr,
                  &m_sluStat,
                  &info,
                  Scalar());
    StatFree(&m_sluStat);

    m_extractedDataAreDirty = true;

    // FIXME how to better check for errors ???
    m_info = info == 0 ? Success : NumericalIssue;
    m_factorizationIsOk = true;
}

template <typename MatrixType> template <typename Rhs, typename Dest> void SuperLU<MatrixType>::_solve_impl(const MatrixBase<Rhs>& b, MatrixBase<Dest>& x) const
{
    eigen_assert(m_factorizationIsOk &&
                 "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");

    const Index rhsCols = b.cols();
    eigen_assert(m_matrix.rows() == b.rows());

    m_sluOptions.Trans = NOTRANS;
    m_sluOptions.Fact = FACTORED;
    m_sluOptions.IterRefine = NOREFINE;

    m_sluFerr.resize(rhsCols);
    m_sluBerr.resize(rhsCols);

    Ref<const Matrix<typename Rhs::Scalar, Dynamic, Dynamic, ColMajor>> b_ref(b);
    Ref<const Matrix<typename Dest::Scalar, Dynamic, Dynamic, ColMajor>> x_ref(x);

    m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
    m_sluX = SluMatrix::Map(x_ref.const_cast_derived());

    typename Rhs::PlainObject b_cpy;
    if (m_sluEqued != 'N')
    {
        b_cpy = b;
        m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
    }

    StatInit(&m_sluStat);
    int info = 0;
    RealScalar recip_pivot_growth, rcond;
    SuperLU_gssvx(&m_sluOptions,
                  &m_sluA,
                  m_q.data(),
                  m_p.data(),
                  &m_sluEtree[0],
                  &m_sluEqued,
                  &m_sluRscale[0],
                  &m_sluCscale[0],
                  &m_sluL,
                  &m_sluU,
                  NULL,
                  0,
                  &m_sluB,
                  &m_sluX,
                  &recip_pivot_growth,
                  &rcond,
                  &m_sluFerr[0],
                  &m_sluBerr[0],
                  &m_sluStat,
                  &info,
                  Scalar());
    StatFree(&m_sluStat);

    if (x.derived().data() != x_ref.data())
        x = x_ref;

    m_info = info == 0 ? Success : NumericalIssue;
}

// the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
//
//  Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
//
//  THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
//  EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
//
template <typename MatrixType, typename Derived> void SuperLUBase<MatrixType, Derived>::extractData() const
{
    eigen_assert(m_factorizationIsOk &&
                 "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
    if (m_extractedDataAreDirty)
    {
        int upper;
        int fsupc, istart, nsupr;
        int lastl = 0, lastu = 0;
        SCformat* Lstore = static_cast<SCformat*>(m_sluL.Store);
        NCformat* Ustore = static_cast<NCformat*>(m_sluU.Store);
        Scalar* SNptr;

        const Index size = m_matrix.rows();
        m_l.resize(size, size);
        m_l.resizeNonZeros(Lstore->nnz);
        m_u.resize(size, size);
        m_u.resizeNonZeros(Ustore->nnz);

        int* Lcol = m_l.outerIndexPtr();
        int* Lrow = m_l.innerIndexPtr();
        Scalar* Lval = m_l.valuePtr();

        int* Ucol = m_u.outerIndexPtr();
        int* Urow = m_u.innerIndexPtr();
        Scalar* Uval = m_u.valuePtr();

        Ucol[0] = 0;
        Ucol[0] = 0;

        /* for each supernode */
        for (int k = 0; k <= Lstore->nsuper; ++k)
        {
            fsupc = L_FST_SUPC(k);
            istart = L_SUB_START(fsupc);
            nsupr = L_SUB_START(fsupc + 1) - istart;
            upper = 1;

            /* for each column in the supernode */
            for (int j = fsupc; j < L_FST_SUPC(k + 1); ++j)
            {
                SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];

                /* Extract U */
                for (int i = U_NZ_START(j); i < U_NZ_START(j + 1); ++i)
                {
                    Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
                    /* Matlab doesn't like explicit zero. */
                    if (Uval[lastu] != 0.0)
                        Urow[lastu++] = U_SUB(i);
                }
                for (int i = 0; i < upper; ++i)
                {
                    /* upper triangle in the supernode */
                    Uval[lastu] = SNptr[i];
                    /* Matlab doesn't like explicit zero. */
                    if (Uval[lastu] != 0.0)
                        Urow[lastu++] = L_SUB(istart + i);
                }
                Ucol[j + 1] = lastu;

                /* Extract L */
                Lval[lastl] = 1.0; /* unit diagonal */
                Lrow[lastl++] = L_SUB(istart + upper - 1);
                for (int i = upper; i < nsupr; ++i)
                {
                    Lval[lastl] = SNptr[i];
                    /* Matlab doesn't like explicit zero. */
                    if (Lval[lastl] != 0.0)
                        Lrow[lastl++] = L_SUB(istart + i);
                }
                Lcol[j + 1] = lastl;

                ++upper;
            } /* for j ... */

        } /* for k ... */

        // squeeze the matrices :
        m_l.resizeNonZeros(lastl);
        m_u.resizeNonZeros(lastu);

        m_extractedDataAreDirty = false;
    }
}

template <typename MatrixType> typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const
{
    eigen_assert(
        m_factorizationIsOk &&
        "The decomposition is not in a valid state for computing the determinant, you must first call either compute() or analyzePattern()/factorize()");

    if (m_extractedDataAreDirty)
        this->extractData();

    Scalar det = Scalar(1);
    for (int j = 0; j < m_u.cols(); ++j)
    {
        if (m_u.outerIndexPtr()[j + 1] - m_u.outerIndexPtr()[j] > 0)
        {
            int lastId = m_u.outerIndexPtr()[j + 1] - 1;
            eigen_assert(m_u.innerIndexPtr()[lastId] <= j);
            if (m_u.innerIndexPtr()[lastId] == j)
                det *= m_u.valuePtr()[lastId];
        }
    }
    if (PermutationMap(m_p.data(), m_p.size()).determinant() * PermutationMap(m_q.data(), m_q.size()).determinant() < 0)
        det = -det;
    if (m_sluEqued != 'N')
        return det / m_sluRscale.prod() / m_sluCscale.prod();
    else
        return det;
}

#ifdef EIGEN_PARSED_BY_DOXYGEN
#define EIGEN_SUPERLU_HAS_ILU
#endif

#ifdef EIGEN_SUPERLU_HAS_ILU

/** \ingroup SuperLUSupport_Module
  * \class SuperILU
  * \brief A sparse direct \b incomplete LU factorization and solver based on the SuperLU library
  *
  * This class allows to solve for an approximate solution of A.X = B sparse linear problems via an incomplete LU factorization
  * using the SuperLU library. This class is aimed to be used as a preconditioner of the iterative linear solvers.
  *
  * \warning This class is only for the 4.x versions of SuperLU. The 3.x and 5.x versions are not supported.
  *
  * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
  *
  * \implsparsesolverconcept
  *
  * \sa \ref TutorialSparseSolverConcept, class IncompleteLUT, class ConjugateGradient, class BiCGSTAB
  */

template <typename _MatrixType> class SuperILU : public SuperLUBase<_MatrixType, SuperILU<_MatrixType>>
{
public:
    typedef SuperLUBase<_MatrixType, SuperILU> Base;
    typedef _MatrixType MatrixType;
    typedef typename Base::Scalar Scalar;
    typedef typename Base::RealScalar RealScalar;

public:
    using Base::_solve_impl;

    SuperILU() : Base() { init(); }

    SuperILU(const MatrixType& matrix) : Base()
    {
        init();
        Base::compute(matrix);
    }

    ~SuperILU() {}

    /** Performs a symbolic decomposition on the sparcity of \a matrix.
      *
      * This function is particularly useful when solving for several problems having the same structure.
      * 
      * \sa factorize()
      */
    void analyzePattern(const MatrixType& matrix) { Base::analyzePattern(matrix); }

    /** Performs a numeric decomposition of \a matrix
      *
      * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
      *
      * \sa analyzePattern()
      */
    void factorize(const MatrixType& matrix);

#ifndef EIGEN_PARSED_BY_DOXYGEN
    /** \internal */
    template <typename Rhs, typename Dest> void _solve_impl(const MatrixBase<Rhs>& b, MatrixBase<Dest>& dest) const;
#endif  // EIGEN_PARSED_BY_DOXYGEN

protected:
    using Base::m_l;
    using Base::m_matrix;
    using Base::m_p;
    using Base::m_q;
    using Base::m_sluA;
    using Base::m_sluB;
    using Base::m_sluBerr;
    using Base::m_sluCscale;
    using Base::m_sluEqued;
    using Base::m_sluEtree;
    using Base::m_sluFerr;
    using Base::m_sluL;
    using Base::m_sluOptions;
    using Base::m_sluRscale;
    using Base::m_sluStat;
    using Base::m_sluU;
    using Base::m_sluX;
    using Base::m_u;

    using Base::m_analysisIsOk;
    using Base::m_extractedDataAreDirty;
    using Base::m_factorizationIsOk;
    using Base::m_info;
    using Base::m_isInitialized;

    void init()
    {
        Base::init();

        ilu_set_default_options(&m_sluOptions);
        m_sluOptions.PrintStat = NO;
        m_sluOptions.ConditionNumber = NO;
        m_sluOptions.Trans = NOTRANS;
        m_sluOptions.ColPerm = MMD_AT_PLUS_A;

        // no attempt to preserve column sum
        m_sluOptions.ILU_MILU = SILU;
        // only basic ILU(k) support -- no direct control over memory consumption
        // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
        // and set ILU_FillFactor to max memory growth
        m_sluOptions.ILU_DropRule = DROP_BASIC;
        m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision() * 10;
    }

private:
    SuperILU(SuperILU&) {}
};

template <typename MatrixType> void SuperILU<MatrixType>::factorize(const MatrixType& a)
{
    eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
    if (!m_analysisIsOk)
    {
        m_info = InvalidInput;
        return;
    }

    this->initFactorization(a);

    int info = 0;
    RealScalar recip_pivot_growth, rcond;

    StatInit(&m_sluStat);
    SuperLU_gsisx(&m_sluOptions,
                  &m_sluA,
                  m_q.data(),
                  m_p.data(),
                  &m_sluEtree[0],
                  &m_sluEqued,
                  &m_sluRscale[0],
                  &m_sluCscale[0],
                  &m_sluL,
                  &m_sluU,
                  NULL,
                  0,
                  &m_sluB,
                  &m_sluX,
                  &recip_pivot_growth,
                  &rcond,
                  &m_sluStat,
                  &info,
                  Scalar());
    StatFree(&m_sluStat);

    // FIXME how to better check for errors ???
    m_info = info == 0 ? Success : NumericalIssue;
    m_factorizationIsOk = true;
}

#ifndef EIGEN_PARSED_BY_DOXYGEN
template <typename MatrixType>
template <typename Rhs, typename Dest>
void SuperILU<MatrixType>::_solve_impl(const MatrixBase<Rhs>& b, MatrixBase<Dest>& x) const
{
    eigen_assert(m_factorizationIsOk &&
                 "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");

    const int rhsCols = b.cols();
    eigen_assert(m_matrix.rows() == b.rows());

    m_sluOptions.Trans = NOTRANS;
    m_sluOptions.Fact = FACTORED;
    m_sluOptions.IterRefine = NOREFINE;

    m_sluFerr.resize(rhsCols);
    m_sluBerr.resize(rhsCols);

    Ref<const Matrix<typename Rhs::Scalar, Dynamic, Dynamic, ColMajor>> b_ref(b);
    Ref<const Matrix<typename Dest::Scalar, Dynamic, Dynamic, ColMajor>> x_ref(x);

    m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
    m_sluX = SluMatrix::Map(x_ref.const_cast_derived());

    typename Rhs::PlainObject b_cpy;
    if (m_sluEqued != 'N')
    {
        b_cpy = b;
        m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
    }

    int info = 0;
    RealScalar recip_pivot_growth, rcond;

    StatInit(&m_sluStat);
    SuperLU_gsisx(&m_sluOptions,
                  &m_sluA,
                  m_q.data(),
                  m_p.data(),
                  &m_sluEtree[0],
                  &m_sluEqued,
                  &m_sluRscale[0],
                  &m_sluCscale[0],
                  &m_sluL,
                  &m_sluU,
                  NULL,
                  0,
                  &m_sluB,
                  &m_sluX,
                  &recip_pivot_growth,
                  &rcond,
                  &m_sluStat,
                  &info,
                  Scalar());
    StatFree(&m_sluStat);

    if (x.derived().data() != x_ref.data())
        x = x_ref;

    m_info = info == 0 ? Success : NumericalIssue;
}
#endif

#endif

}  // end namespace Eigen

#endif  // EIGEN_SUPERLUSUPPORT_H
